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Abstract. In this paper we study the reconstruction of a bandlimited signal 
from samples generated by the integrate and fire model. This sampler allows 
us to trade complexity in the reconstruction algorithms for simple hardware 
implementations, and is specially convenient in situations where the sampling 
device is limited in terms of power, area and bandwidth. 

Although perfect reconstruction for this sampler is impossible, we give a 
general approximate reconstruction procedure and bound the corresponding 
error. We also show the performance of the proposed algorithm through nu- 
merical simulations. 

Keywords: integrate and fire, non-uniform sampling, bandlimited func- 
tion. 



1. Introduction 

The integrate and fire (IF) model is well known in computational neuroscience 
as a simplified model of a neuron [8j [12] and is typically used to study the dynam- 
ics of large populations. The model consists of a leaky integrator followed by a 
comparator. The leak corresponds to a gradual loss of the value of the integral. 

More recently, the IF model has also been considered as a sampler [4, 16, 10] ITT], 
where the sampler output is tuned to the variation of the integral of the signal. 
This feature can be exploited when sampling neural recordings, for which relevant 
information is localized in small intervals where the signal has a high amplitude [3] . 

The block diagram of the sampler is presented in Figure [T] At every instant s, 
the continuous input x(t) is integrated against an averaging function Uk, s (t) and 
the result is compared to a positive and negative threshold. When either of these is 
reached, a pulse is created at time tk = s representing the threshold value (positive 
or negative), the value of the integrator is then reset and the process repeats. The 
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Figure 1. Block diagram for the BIF model. 



output is a nonuniformly spaced pulse train, where each of the pulses is either 
1 or -1. The averaging function Uk, s (t) is defined by e a ^~ 3 ^X[ tki3 ^ where Xj is 
the characteristic function of I and a > is a constant that models the leakage 
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of the integrator due to practical implementations. The precise firing condition 
determining the pulses is: 



The simplicity of the sampler translates into an efficient hardware implementation 
which saves both power and area when compared to conventional analog-to-digital 
converters (ADC) [4J. These constraints are severe, in the case of wireless brain 
machine interfaces [13 j, for which the entire system has to be embedded inside the 
subject. Hence, the IF sampler allows us to move the complexity of the design into 
the reconstruction algorithm while providing a simple front end at the sampling 
stage. 

The problem of reconstructing a signal from the IF output should be distin- 
guished from the study of the dynamics of a population of neurons, when some 
stochastic assumption is made on the firing parameters [9] . In this article we study 
the deterministic reconstruction of a bandlimited signal from the integrate and fire 
output. Part of the challenge of this stems from the fact that the sampling map that 
associates a function to its samples is non-linear. Indeed, we see from Equation (pQ) 
that the magnitude of the samples is always 6. Moreover, exact reconstruction for 
the IF sampler is impossible since the output of the sampler does not completely 
determine the signal (see Example Q] below.) 

In this article, we will show that it is however possible to approximately recon- 
struct a bandlimited signal in L°° norm with an error comparable to the threshold 6. 
Moreover, we give a concrete reconstruction procedure which is of course non-linear 
but, nevertheless, easy to implement. Since in many situations the IF sampler is so 
much more convenient to implement than conventional analog-to-digital converters, 
the loss of accuracy in the reconstruction is a very reasonable trade-off [4 , specially 
if the final analysis of the reconstructed data tolerates some small error [3]. 

The methods considered so far [11] reconstruct the signal / from the system 
of equations (f,Uk) — ±0 (cf. Equation (pQ)), thus treating the reconstruction as 
a (linear) average sampling problem (see [3 [TJ [151 El ) • These approaches impose 
density restrictions on the set of sampling functions {uk} k (cf. Equation (pQ).) Since 
these sampling functions depend on the signal, the density constraints on them are 
somehow unnatural. 

The key for the reconstruction method that we develop lies in the observation 
that the information derived from the IF output is much richer than the mere 
system of equations (/, i/fc) = ±0. It also contains the information that no proper 
subinterval [tk,t'] of tfc+i] satisfies Equation (pp). We will exploit this extra infor- 
mation to give an approximate reconstruction procedure for a general bandlimited 
function. Since the sampling process starts at a certain instant to, an additional 
assumption on the size of / before to is required in order to fully reconstruct /. 
Roughly speaking, the assumption means that the sampling scheme would not have 
produced any pulse before to. 

In Section [2] we formally describe the output of the IF sampling scheme. This 
output depends on an initial time to when the process is started and two parameters: 
the threshold and the constant a > modeling the leakage on the sampler. In 
Section [3] we show that the IF output is always a finite sequence. Section [4] gives 
the approximate reconstruction procedure and Section [5] presents some numerical 
experiments. 
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2. The integrate and fire sampling problem 

We now define precisely the integrate and fire sampling scheme. Throughout the 
article we will assume the following. 

Assumption 1. A bandlimited function f G PWq and numbers to G R ; a, > 
are given. 

Here, PWq is the Paley- Wiener space 

PWn :={fe L 2 (R) | supp(/) C [-0, Q] }, 

of (complex- valued) bandlimited functions and f(w) := J R f(x)e~ 2n,lwx dx is the 
Fourier transform of /. We call to the initial time, a the firing parameter and the 
threshold. Using these parameters we formally define the output of the sampler. 
We first define recursively a finite or countable sequence to < ... <tj . .. called the 
time instants. Suppose that the instants to < . . . < tj have already been defined 
and consider the function Fj : [tj, +oo) — > C given by 



Observe that Fj is continuous and Fj(tj) = 0. If \Fj(t)\ < 0, for all t > tj, then the 
process stops. If \Fj(t)\ > 0, for some t >tj, by the continuity of Fj, we can define 
tj+i as the minimum number satisfying the equation 



(2) 



f(x)e aix ~ t ^ l) dx 



Clearly, in this case tj+i > tj. 

We have defined a finite or countable sequence of points to < ... < tj — We 
will prove in Proposition [T] that this sequence is in fact finite. Let us assume the 
time instants {to, • • • ,t n } and define the samples {qi, . . . , q n } by, 

(3) qj := f(x)e^ x -^dx, (1 < j < n). 



Observe that, by the definition of the time intervals, \qj\ = 0. 

The output of the sampler is formally given by the time instants {to, . . . ,t n } 
and the numbers {qi, . . . ,q n }. We say that this output has been produced by 
the integrate and fire. The succeeding results apply generally to complex- valued 
functions, but in the case of the application that motivated this sampling scheme, 
the signal is taken to be real-valued, and the output of the sampler is encoded as a 
train of impulses, where only the sign of the samples qj is stored. 

3. Some remarks on the IF output 

First we note that bandlimited functions are not completely determined by the 
output of the IF sampler. 

Example 1. There are non-zero bandlimited signals that will never produce an 
output from the sampler. Take for instance fo(x) = ^2~^ • Since fo(u) = 



6 
2 

rt 



max{l — \co\, 0} ; fe is bandlimited. We have for any to Gl, J t fe(x)e a ^ x ^dx 



< 



SL dx<lf R s -^dx=l<6,forallt>t . 



sin (n x) 
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We now prove that the set of time instants produced by the IF sampler is indeed 
finite and give some bounds on its distribution. To this end we introduce some 
auxiliary functions that will be used throughout the remainder of the article. 

Consider the function g : R — > R given by g(x) = e~ ax X[o y00 ] an d define, 

(4) v{t) := (/ * g)(t) := /* f(x)e a ^dx. 



Since g G ^(R), v G PWq. In the Fourier domain, v and / are related by 

(5) f(w) = (2i\iw + a) v(w). 
In the time domain, this can be expressed as 

(6) m = ^ + av{t ). 

Observation 1. The function v is continuous and v(t) — > 0, when t — ► ±oo. 

Proof. We have already observed that v G PWq. Since v G L 2 and supp(v) C 
[—17, fi], we have that v G L 1 and the conclusion follows from the Riemann-Lebesgue 
Lemma. □ 

The following straightforward equation relates v to the integrate and fire process. 

(7) / f(x)e a{x - t) dx = v(t) - e^-^vis), s<t. 

J s 

We can now prove that the output of the IF process is finite. 

Proposition 1. Under Assumption 1, the following holds. 

(a) The set of time instants produced by the integrate and fire scheme is a finite 
set {t , . . . ,t n }. 

(b) The numbers of time instants tj in a given finite interval [a, b] is bounded 
by 

tt (6 _ a) l/2 + L 



(c) If f is integrable, the total number of time instants is bounded by 

ii/iii 

' 

Proof. We first prove (b) and (c). Let [a, b] be an interval and let {tj : . . . , t J+m _i} 
be m consecutive time instants contained in [a, b]. If m < 1 the bound is trivial, so 
assume that m > 2. For each < k < m — 2, using Equation ([2j) we have, 

( 3+k+1 f(x)e a{x ~ t ^ k + l) dx 

Z 3 + k 

/tj + k + 1 
\m\dx. 
-j + k 

Summing over the m — 1 intervals determined by the points {tj, . . . ,i/ +m _i} we 
have, 

(8) (ra-l)<9< / \f(x)\dx. 

J a 
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Letting a — — oo and b = +00 yields (b). For (a), Holder's inequality gives, 

(m-l)6< \\fh{b-af/\ 

and the conclusion follows. 

Now we prove (a). Assume on the contrary that the IF process goes on forever 
producing an infinite set of instants {tj : j > 0}. Given s > to, by part (b), only 
a finite number of instants tj belong to [to,s]. Therefore t n — > +00, as n — > +00. 
Using Equations (j7j) and ([2]) it follows that, 

^ f(x)e a{x - t ^ l) dx 

= v(t j+1 )-e a ^-^Mtj) <K*i+i)l + Kti)|. 
This contradicts Observation [U □ 



4. The reconstruction 

We now address the problem of approximately reconstructing a bandlimited 
function from the integrate and fire output. Since the samples are taken in the 
half-line [to, +00) we will make some assumption about the size of / before the 
initial instant. Roughly speaking, the integrate and fire process would not have 
produced any sample in the interval (— 00, to]. 

Assumption 2. The function defined in Equation @ satisfies, 

\v(t)\ < 0, for allt < t . 

Note that by Observation [TJ any to <C satisfies this assumption. To approxi- 
mately reconstruct / we will first approximately reconstruct v from the integrate 
and fire output and then derive information about / by means of Equation ([5j). 
We will use the structure of the IF process to produce a number of approximate 
samples for v. 

First we argue that, from the output of the IF process, we have enough infor- 
mation to approximate v on the time instants {to, . . . , t n }. Rewriting Equation ([3]) 
in terms of v (cf. Equation Q) we have, 

(9) v(t j+1 )=e a ^- t ^v(t j ) + q j+1 , (0<j<n-l). 

Since the value v(to) may not be exactly known we cannot determine from this re- 
currence relation all the values v(tj). However, we can construct an approximation 
to these values. Let wq := and define recursively, 

(10) w j+1 = e^-^Wj + q j+u (0<j<n- 1). 

Observe that Assumption [2] implies that \wo — v(to)\ < 0. Using this estimate as a 
starting point we can iterate on Equation ([9|) and ([10]) to get, 

(11) \wj-v(tj)\ <0, (0< j<n). 

Consequently, using only the output of the IF sampling scheme, we have constructed 
a set of values {i^o, . . . ,w n } that approximates v on the instants {to, . . . , t n }. The 
second step is to approximate v on an arbitrary point of R. 
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To this end observe that, according to the definition of tj as the minimum number 
satisfying Equation ([2j), we have that, 



f(x)e a ^ x - t) dx 



< 0, for all t E [tj,tj+i\. 



Rewriting this inequality in terms of v (cf. Equation 0) gives, 



(12) 



v(t) 



< 9, for all t G [tj,tj+i]. 



Combining this last inequality with (fTT|) yields, 



(13) 



v(t)-e a ^-^w 



< 2(9, for all t G [tj,t j+1 ]. 

We now show that this inequality allows us to approximate v anywhere on the line. 

Claim 1. Given an arbitrary time instant t G M, choose x G M d in the following 
way: 

(a) if t < to, let x := ; 

(b) if t belongs to some (unique) interval [tj,tj+i), let x := e a ^ tj ~^Wj , 

(c) ift> t n , let x := e a ( tn ~^w n . 
Then, \v(t) - x\ < 20. 

Remark 1. Observe that the procedure to obtain x from t depends only on the 
output of the IF process. 

Proof. For case (a), the conclusion follows from Assumption [2l For case (b), the 
conclusion follows from Inequality ([T3j) . For case (c), the fact that the fire condition 
is never satisfied after t n gives, 



(14) 



«(*) 



*(t n -t) 



v(t n ) 



< 



□ 



Combining this estimate with Inequality (fTT|) . the conclusion follows. 
We will now choose a window function. 

Assumption 3. A Schwartz class function ip such that 

• i[j = 1 on [— 0,0] ; and, 

• i/j is compactly supported, 
has been chosen. 

Since v G PWq, the classic oversampling trick for bandlimited functions (see for 
example [6] or [7]) implies that there exists a number < /3 < (20) _1 such that 



(15) 



v = ^2v(j3k)il>{--l3k). 



Using the procedure described in Claim HJ we produce a set {sk} ke % such that 

(16) \v(f3k) -s k \< 2(9, for all k G Z. 
Let (f be the function defined by, 

(17) — i^iriw + a) ^(w). 
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It follows that (p is also a Schwartz function. Moreover, using Equation ([5]) we have 
that, 

(18) f = J2v(PkM--(3k). 

kez 

Observe that, since v G PWq, the sequence {v(/3k)} k G £ 2 and the series in Equa- 
tion ([18)) converges in L 2 and uniformly - in fact, it converges in the Wiener amal- 
gam norm W(Co, L 2 ), see for example [6], [7] and [2J.) 

Now we can define the approximation of / constructed from the IF samples. Let, 

(19) 

Since, by Inequality ([16]) . the sequence {sk} k is bounded and ip is a Schwartz 
function, it follows that Equation ([19]) defines a bounded function and that the 
convergence is uniform (see [6] or [2].) 

The reconstruction algorithm consists then of calculating the approximated sam- 
ples {sk} k following Claim 1 and then convolving them with the kernel (/?, that can 
be pre-calculated. 

We now give a precise error bound for the reconstruction. 

Theorem 1. Under Assumptions 1, 2 and 3, the function defined by Equation ([19]) 
satisfies, 

11/ - /Hoc < co, 

for some constant C that only depends on Q and the window function chosen in 
Assumption^ 

Proof. According to Equations ([18]) . ([19]) and Inequality ([16]) . 

11/ - /Hoc < supess^ \v(J3k) - s k \ \ip{. - (3k)\ 
kez 

< 2<9supess^ \(p(- -0k)\. 

kez 

It suffices to define C := 2sup^ fcGZ \cp(- — /3k)\. Since ip is a Schwartz function, 
C < +oo (see for example [5]). □ 

Remark 2. We currently do not know what choice of the window function ip min- 
imizes the constant in the theorem. A more detailed study of the choice of the 
window function should not only consider the size of that constant but also the rate 
of convergence of the series in Equation ([19]) . 

5. Numerical experiments 

We study the behavior of the reconstruction algorithm under variations in the 
threshold and the oversampling period for a specific choice of reconstruction ker- 
nel ip. The test signal / is of finite length and real valued, produced as a linear 
combination of five 'sine' kernels (sin(7nr)/(7nc)) at a 1Hz frequency, with random 
locations and weights. The amplitude of the input has been normalized to 1. Al- 
though the theory covers infinite dimensional spaces, our simulations are limited 
by the practical implementations of the sampler and algorithms. The effects of 
truncation and quantization are not considered here. 
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This signal is encoded by the IF sampler with a = 1 and recovered using the 
procedure described in Section [U The reconstruction kernel ^ is a raised cosine, 
defined by, 

(20) = smc(t/T s ) cos(njt/T s ) (l - 

where 7 = 0.5 and T s = 0.25 are determined by the maximum input frequency 
Q and the desired oversampling period /3 (cf. Equation (fT5|) .) Figure 2(a)| shows 



the raised cosine x/j in the time and frequency domain. Observe that the spectrum 
of ip is constant for frequencies less than the input bandwidth and then decays 
smoothly towards zero. The corresponding kernel <p (cf. Equation ([T7|) ) is shown 
in Figure [2(b)] Using ip we recover / (cf. Equation ([19]) .). an approximation of / 




Figure 2. Reconstruction kernels. 



as shown in Figure [3j As expected the error decreases in regions with high density 
of samples. This behavior is evident from Figure [H the dense regions imply that 
the uniform samples will most likely coincide with the estimated values of v(t) 
at the sample locations. On the other hand, for samples that are far apart the 
approximation follows a exponential decay from its original value which is not the 
natural trend in the signal. Figure [4] shows v(t) (solid line), and the approximated 
samples of v on the lattice /?Z, constructed using the procedure described in Claim 
Q] called {sk} k and the envelope v(t) ± where these samples are known to lie 
(dashed line.) 

Currently the reconstruction algorithm uses the approximated samples of v(t) 
at the pulse locations to define the piecewise exponential bound and estimate the 
reconstruction coefficients on the uniform lattice. Based on the numerical experi- 
ments the algorithm can be improved by including the estimated value of v(t) at 
the pulse locations although it implies reconstruction on a nonuniform grid. For 
both cases similar error bounds can be defined as in Theorem [TJ The variation of 
the error in relation to the threshold (pulse rate) is shown in Figure The error 
depends on the choice of generator and the oversampling period /?, as seen in Fig- 
ure [U The relationship between the kernels and the optimal oversampling period 
is still not evident. 
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Figure 3. Reconstruction of f(t) from the impulse train. 
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Figure 4. Reconstruction of v(t) with (3 = 1/4, 6 = 0.05 
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Figure 5. Variation of the error ||/ — /||oo i n relation to the 
threshold and pulse rate (dotted line) with /3 = 0.25. 




Figure 6. Variation of the error in relation to the oversampling 
period for different thresholds. The error is defined as ||/ — /||oo- 
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